namespace GAUSS {
	constexpr double eps = 1e-9;

	// 一个n * (n + 1)的矩阵，最后一列是方程的值
	// 求方程解
	vector<double> gauss(vector<vector<double>> a) {
		if (!a.size()) {
			return {};
		}
		assert(a.size() + 1 == a[0].size());

		int n = (int)a.size(), r = 0;
		for (int c = 0; c < n && c < n; c++) {
			int t = r;
			for (int i = r + 1; i < n; i++) {
				if (fabs(a[i][c]) > fabs(a[t][c])) {
					t = i;
				}
			}
			if (fabs(a[t][c]) < eps) {
				continue;
			}
			for (int i = c; i <= n; i++) {
				swap(a[t][i], a[r][i]);
			}
			for (int i = n; i >= c; i--) {
				a[r][i] /= a[r][c];
			}
			for (int i = r + 1; i < n; i++) {
				if (fabs(a[i][c]) < eps) {
					continue;
				}
				for (int j = n; j >= c; j--) {
					a[i][j] -= a[r][j] * a[i][c];
				}
			}
			r++;
		}
		
		if (r < n) {
			for (int i = r; i < n; i++) {
				if (fabs(a[i][n]) > eps) {
					// 无解
					return {-1}; 
				}
			}
			// 不是唯一解
			return {};
		}
		for (int i = n - 2; i >= 0; i--) {
			for (int j = i + 1; j < n; j++) {
				a[i][n] -= a[j][n] * a[i][j];
			}
		}
		vector<double> ret;
		for (int i = 0; i < n; i++) {
			ret.push_back(a[i][n]);
		}
		return ret;
	}

	// 求行列式值
	int gauss(vector<vector<int>> a, int p) {
		int ret = 1, n = a.size();
		for (int i = 0; i < n; i++) {
			for (int j = i + 1; j < n; j++) {
				while (a[i][i]) {
					int x = a[j][i] / a[i][i];
					for (int k = i; k < n; k++) {
						a[j][k] = (a[j][k] - 1ll * x * a[i][k] % p) % p;
						a[j][k] = (a[j][k] + p) % p;
					}
					a[i].swap(a[j]);
					ret = -ret;
				}
				a[i].swap(a[j]);
				ret = -ret;
			}
		}
		for (int i = 0; i < n; i++) {
			ret = 1ll * ret * a[i][i] % p;
		}
		return (ret + p) % p;
	}
}